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Abstract 

Background: Variations in gene expression, mediated by epigenetic mechanisms, may cause broad phenotypic 
effects in animals. However, it has been debated to what extent expression variation and epigenetic modifications, 
such as patterns of DNA methylation, are transferred across generations, and therefore it is uncertain what role 
epigenetic variation may play in adaptation. 

Results: In Red Junglefowl, ancestor of domestic chickens, gene expression and methylation profiles in thalamus/ 
hypothalamus differed substantially from that of a domesticated egg laying breed. Expression as well as 
methylation differences were largely maintained in the offspring, demonstrating reliable inheritance of epigenetic 
variation. Some of the inherited methylation differences were tissue-specific, and the differential methylation at 
specific loci were little changed after eight generations of intercrossing between Red Junglefowl and domesticated 
laying hens. There was an over-representation of differentially expressed and methylated genes in selective sweep 
regions associated with chicken domestication. 

Conclusions: Our results show that epigenetic variation is inherited in chickens, and we suggest that selection of 
favourable epigenomes, either by selection of genotypes affecting epigenetic states, or by selection of methylation 
states which are inherited independently of sequence differences, may have been an important aspect of chicken 
domestication. 
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Background 

Chickens were domesticated from the Red Junglefowl 
(RJF) about 8000 years ago [1,2], and the changes in mor- 
phology, physiology and behaviour as a response to this 
have been immense. For example, most domesticated 
chickens grow to at least twice the size of RJF, become 
sexually mature at a lower age, lay manifold more and lar- 
ger eggs, show a wide variation in plumage colour and 
structure, and have a different behaviour in a number of 
contexts, such as reduced fearfulness [3-6]. In general, 
domestic animals are assumed to have adapted to a life 
among humans by evolving higher flexibility in diet, better 
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ability to breed in captivity, less stress susceptibility, and a 
more socially tolerant disposition [5,6]. It has been sug- 
gested that epigenetic mechanisms might be involved in 
cases like this [7] where wide-encompassing phenotypic 
changes occur in a short evolutionary time. 

However, there is limited knowledge of the extent to 
which expression and epigenetic profiles are inherited in 
animals. Reliable inheritance is necessary in order for 
epigenetic variation to be a major component of any evo- 
lutionary process. We have earlier shown that stress- 
induced modifications in both behaviour and brain gene 
expression profiles in domestic chickens are to some 
extent transferred to the offspring [8,9], and other studies 
have shown similar transgenerational transmission in 
other species, including humans [10-12]. This indicates 
that some epigenetic variation may indeed be inherited, 
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but the details and significance of this, as well as its puta- 
tive evolutionary significance, remain to be elucidated. 

One of the possible epigenetic mechanisms, which 
could be related to variation in gene expression, is 
methylation of cytosine, preferentially in so called CpG- 
islands of promoter regions [13,14]. Therefore, we 
targeted methylation and gene expression simulta- 
neously to investigate whether any of those, or both, 
would differ between two populations of chickens, 
recently separated by domestication. We hypothesised 
that both methylation and gene expression would differ 
between the populations and show transgenerational 
stability, opening the possibility for both to be involved 
in domestication-related phenotypic changes. 

By using expression and methylation arrays on hypotha- 
lamus samples, we show that profiles of gene expression 
as well as promoter methylation differ between domesti- 
cated White Leghorn layer chickens and their ancestors, 
the Red Junglefowl. There were also similar differences, 
although less pronounced, between phenotypically differ- 
ent families within breeds. The differences were largely 
maintained in the offspring, demonstrating a reliable 
inheritance of epigenetic states, and for some of the genes 
the differential methylation was maintained after eight 
generations of intercross. Our results therefore suggest 
that selection of favourable epigenetic variants may have 
been an important aspect of chicken domestication. 

Results and discussion 

Brain gene expression differences within and between 
populations 

In this experiment, we studied variations in gene expres- 
sion and methylation in brains of RJF and domesticated 
White Leghorns (WL), and their offspring. We focussed 
on thalamus and hypothalamus, brain regions involved in 
fear and stress responses, both of which have changed 
significantly during domestication [3,6]. Within each 
population, we selected parental animals with divergent 
phenotypes in order to maximise the within population 
genetic variation. Specifically, we used two pairs of each 
population, with pairs within population differing in their 
behaviour in a series of previously validated tests of stress 
reactions in chickens [6,15]. From these, totally 73 off- 
spring were hatched and reared until three weeks of age, 
when they were tested in a fear test, similar to that used 
in the parents. 

In both breeds, body weights differed between families 
in both generations, and behavioural scores, as measured 
in the fear tests, differed between families in both genera- 
tions of WL, but not RJF (Additional file 1). Hence, mor- 
phological, and to some extent behavioural, phenotypes 
showed a significant and transgenerationally stable varia- 
tion in the animals used for the present study. It should 
be noted that phenotyping was done at different ages in 



the two generations, which may have been the reason for 
the lack of transgenerational correlation in fear behaviour 
in RJF. All eight parents were sacrificed at an age of 373 
days, and 48 offspring (12 from each pair) at 21 days, and 
from each brain, the thalamus-hypothalamus region was 
removed for extraction of both DNA and mRNA. For the 
offspring, eight pools of both were prepared, each con- 
sisting of six same-sex samples within families. Hence, 
there were in total eight parental single-animal samples, 
and eight pools of offspring samples. The mRNA was 
hybridized to a 38K Affymetrix chicken gene expression 
microarray, and the DNA was used for subsequent tiling 
array analysis of methylation. Between populations, there 
were in total 281 significantly (FDR-corrected P < 0.05) 
differentially expressed (DE) genes in the parents, and 
1674 in the offspring. The lower number of DE genes in 
the parents could possibly be an effect of the lower 
power of detection given the smaller biological sample 
size in this generation. Between families within popula- 
tions, only a few genes were significantly DE, and also 
DM was less frequent between families (Additional file 
2). This indicates that expression and methylation pro- 
files are relatively stable within breeds, but both may 
have changed considerably during domestication. 

Transgenerational stability of gene expression profiles 

Out of the significantly DE genes in the parents (compar- 
ing populations), 86% percent (n = 242) were also signifi- 
cantly DE in the offspring (Additional file 3), and there 
was a distinct similarity in the expression differences in 
both generations (Figure la). The overall pattern of fold- 
change levels between populations (regardless of whether 
they were significant) was strongly correlated over genera- 
tions (Figure 2 a), further showing a transgenerational sta- 
bility in gene expression profiles. Also within populations, 
the overall pattern of fold-change levels between families 
was highly correlated across generations (Figure 2 b-c). 
We further used signalling intensities of individual pro- 
besets on each microarray to correlate global expression 
levels between parents and their own offspring, compared 
to offspring of other birds, and found a significantly 
higher correlation within families than between (mean 
difference in correlation coefficients 0.0017 ± 0.0002 
(SEM), t = 8.2, P < 0.001). This was true both for RJF and 
WL, and further supports that specific brain gene expres- 
sion profiles are indeed inherited. 

Gene methylation: inheritance and differences between 
populations 

For analysis of differential methylation (DM), we 
selected 3623 genes from the list of genes which had the 
highest fold changes in DE in both generations, both in 
the between- and within-population comparison. Note 
that only 281 of these were significant in parents and 
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Figure 1 Gene expression and methylation differences between populations and across generations, a. Heat map, showing the 
clustering of 242 differentially expressed genes, comparing parental Red Junglefowl (RJF) and White Leghorn (WL) layers, and their offspring, b. 
Heat map, showing clusters of differentially methylated genes comparing parental RJF and WL, and their offspring. Note that the gene set in b is 
not the same as in a. c. Average methylation levels (signal intensity from the microarray) for RJF and WL in the 145 promoter regions included 
in the heatmap in panel b. Circle-diagram displays the percentage of the promoter regions which were hypermethylated (green) and 
hypomethylated (red) in White Leghorns. 



1674 in offspring. For each of these genes, 50-75 bp- 
probes representing a region spanning from -7.25 kb 
upstream to +3.25 kb downstream of the transcription 
start point (hence mostly covering promoter regions and 
other cis-acting regulatory elements) were placed on a 
custom made tiling array. Methylated DNA immune 
precipitation (MeDIP) was used to enrich methylated 
DNA fragments, and after labelling and hybridisation, 
the relative levels of methylated to un-methylated DNA 
was assessed for each probe. 

Out of the 3623 selected genes, 239 were significantly 
DM (FDR-corrected P < 0.05) when comparing RJF and 
WL parents, and 821 were DM in the corresponding 
comparison in the offspring. A smaller number were 



classified as DM when comparing between families 
within population (Table S2). A heat map of the genes 
classified as DM in both generations showed a highly 
consistent pattern across generations (Figure lb). 
Furthermore, DM levels were significantly correlated 
between generations when comparing RJF with WL 
(Figure 2 d), and also to a lesser degree when comparing 
WL, but not RJF families (Figure 2 e-f). 

Of the 145 genes which were significantly DM in both 
generations (Additional file 4; Additional file 5), 79% 
were hypermethylated in WL (Figure lc). This is a 
highly significant bias (% 2 = 49.8, P < 0.0001), indicating 
that this breed has acquired novel methylation patterns 
during its selection history. 
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Figure 2 Correlations of differential expression and methylation of genes between generations, a-c, Correlations between generations of 

differential gene expression, comparing Red Junglefowl and White Leghorn (a); families within Red Junglefowl (b); and families within White 

Leghorns (c); d-f, Correlation between generations of differential methylation, comparing Red Junglefowl and White Leghorn (d); families within 

Red Junglefowl (e); and families within White Leghorns (f). The genes included in the analyses were selected based on the fact that they 

simultaneously occurred in both parents and offspring among top 1000 DE or DM, sorted by fold change. 
\ J 



We further analysed the relationship between DM and 
DE on the 3623 selected genes. There was no overall cor- 
relation between the level of DM of a gene (% of DM 
probes) and the degree of DE of the same gene (Additional 
file 4). Furthermore, there was no overrepresentation of 
DE genes among the top 100 DM promoters when com- 
pared to a random sample of 100 DM genes (% 2 = 2.1, P > 
0.05). This is contrary to the common notion that methy- 
lation causes down-regulation of gene expression, but 
similar findings have recently been reported from other 
species, for example humans [16,17]. The finding is quite 
surprising, and indicates that the specific sites of 



methylation may be of major importance for gene regula- 
tion. For example, there may be a substantial difference 
between methylation of transcription factors compared to 
insulator sequences. Since we only analysed a 10 kb region 
around the transcription start site of each gene, we can 
not exclude that DM in other, more distant regulatory 
regions may be more closely connected to the expression 
level. 

To illustrate examples of the transgenerationally 
stable methylation patterns observed, we show methy- 
lation graphs for four genes (ABHD7, GAB1, KSR1 and 
PCDHAC1) in Figure 3. In all four, the methylation 
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Figure 3 Transgenerational stability of methylation patterns in specific genes a-h, Differential methylation levels (Log2 fold change) of 
promoter regions, comparing Red Junglefowl and White Leghorn, are shown with a resolution of 50-75 bp-regions in parents and offspring for 
each of the genes a, b, ABHD7, c, d, GABl, e, f, KSRl, and g, h, RAPGEF1. Transcription direction and exons (blue boxes) are shown, and red 
arrows point at bars or groups of bars where significant levels of differential methylation are found (also indicated by red bars). Locations of 
CpG-islands are indicated in yellow for each region. 
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pattern was reliably inherited, shown by the fact that 
the DM pattern was highly similar in parents and off- 
spring. ABHD7 showed extensive DM ranging several 
kb downstream of the transcription start site. In none 
of the four genes, the significantly DM loci were in 
CpG-islands, so methylation must have targeted cyto- 
sines in other genomic contexts. Extensive methylation 
of non-CpG regions have recently also been reported 
for the human methylome [16,17], and it remains 
unknown which functions these epigenetic variants 
may serve. 

Verification of differential methylation with independent 
animals, tissues and method 

To verify the results of the array-based methylation analy- 
sis, we arbitrarily selected four genes, which were DM on 
the tiling arrays in either parents or offspring, FUCA1, 
PCDHAC1, TXNDC16, and RUFY3, and replicated the 
findings for those, using a different technique and a differ- 
ent animal material. Hypothalamus /thalamus regions from 
eight five-weeks old RJF and eight WL (same strains as 
earlier, but different parents) were dissected and treated as 
described above. The DNA was bisulfite-treated, and the 
degree of methylation was determined in the regions that 
were significantly DM on the tiling array using methyla- 
tion sensitive high resolution melting (MS-HRM) analysis. 

All four genes were significantly DM in the same direc- 
tion as found on the tiling array (FUCA1 and PCDHAC1 
hypermethylated in WL; RUFY3 and TXNDC16 hypo- 
methylated) (Figure 4a). This suggests that the tiling array 
produced reliable results and that the observed methyla- 
tion differences are representative for the population dif- 
ferences at large. 

In order to check for tissue-specificity of the DM, we 
also performed HRM on the same four genes, using 
DNA-pools prepared from cerebellum and blood from 
the offspring samples included in the tiling arrays. All 
four genes were significantly DM in cerebellum. In blood, 
FUCA1 and PCDHAC1 were significant, while RUFY3 
showed a tendency for DM (P = 0.08) (Figure 4 a). The 
fact that TXNDC16 was not DM in blood indicates that 
this gene shows tissue-specific, heritable methylation. 

Genetic stability of methylation differences 

There is a risk that the methylation differences detected 
by the MeDIP technique could be affected by sequence 
differences in the promoter regions used for the arrays. 
To exclude this possibility, we used the recently pub- 
lished resequencing data of Red Junglefowl and domes- 
tic chickens [18] to check the 145 significantly DM 
probes in both parents and offspring for possible dele- 
tions, insertions and SNP density. Apart from occasional 
SNPs (Additional file 5), no major sequence differences 
were detected. 



The methylation differences observed may be a result 
of either inheritance of the epigenetic changes indepen- 
dently of genetic changes, or result from sequence dif- 
ferences which secondarily affect methylation at close or 
remote loci. This is more difficult to differentiate, since 
it would require extensive resequencing data of the indi- 
viduals actually used in the study, combined with, for 
example, methylation QTL-studies. 

To suggestively analyse whether differential methylation 
of specific loci are caused by sequence differences we 
decided to study its genetic stability and segregation over 
several generations. For this purpose, we used a total of 18 
birds from the eighth generation of an intercross between 
RJF and WL. In this population, genetic recombinations in 
each generation have broken up the linkage between adja- 
cent loci, and we could therefore check for both stability 
of the methylation sites, and for possible cis- or trans-reg- 
ulation of these. 

From this group of advanced intercross birds, we 
selected individuals, which were homozygous for either 
the WL or RJF-allele, or heterozygotes, of an SNP located 
within 176-1449 kb of the locus showing DM. Using HRM 
analysis on DNA from blood, extracted from these differ- 
ent genotypes, we again analysed the methylation on 
FUCA1, PCDHAC1 and RUFY3 in these individuals (Fig- 
ure 4b). For FUCA1, we found two different non-signifi- 
cant, but distinct, methylation levels, where the birds 
homozygous for the WL-marker were hypermethylated, 
and heterozygotes were similar to the ones homozygous 
for the RJF-marker (P = 0.07). With respect to PCDHAC1, 
the three genotypes were significantly different (P < 0.001), 
with heterozygotes having a methylation level falling 
between the hypermethylated WL homozygotes, and the 
RJF genotypes. RUFY3 showed a high level of methylation, 
which was not significantly different between the three 
genotypes. Hence, two of the three DM loci were stable 
over the eight generations of intercrossing, and tended to 
segregate according to genotype at the locus. This is con- 
sistent with a cis-regulating mechanism, showing a domi- 
nant inheritance of hypomethylation in genotypes with 
RJF alleles for FUCA1, and an intermediate, codominant 
inheritance in PCDHAC1. RUFY3 may possibly be under 
control of trans-acting loci, which have segregated during 
the intercrossing. 

Although these results are not conclusive, they suggest 
that sequence differences may determine the DM for at 
least two of the three loci, possibly for all of them. This 
further suggests that selection during domestication may 
have targeted genotypes which modify the epigenomes, 
perhaps affecting phenotypes indirectly. 

Genetic pathways 

To examine which genetic pathways and functions that 
may have been affected by DE and DM, we performed a 
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Gene Tissue 



RJF 
MeaniS.E. 



WL 
Meant S.E. 



df* 



SNP # 



FUCA1 



Hypothalmus/Thalamus 77.35 8 26.9910.34 8 43.6314.32 -3.84 7.09 <0.01 

Cerebellum 77.4 7 23.3910.68 7 38.87+4.57 -3.38 6.04 <0.05 

Blood 77.25 6 17.85+0.24 6 31.2513.61 3.71 5.04 <0.05 



No 



PCDHAC1 

Hypothalmus/Thalamus 79.3 

Cerebellum 79.1 

Blood 79 



9.17+ 1.13 
9.1611.05 
8.69+ 1.96 



25.171 1.34 
20.82+ 1.35 
16.6910.80 



-9.13 
-6.83 
3.2 



13.61 O.OOOOl 
11.31 <0.0001 
6.13 <0.05 



No 



RUFY3 

Hypothalmus/Thalamus 76.05 

Cerebellum 76.05 

Blood 75.9 

TXNDC16 

Hypothalmus/Thalamus 75.35 

Cerebellum 75.35 

Blood 75.2 



57.71+0.39 
57.4610.26 
53.42+0.43 



46.78+0.76 
50.61+0.46 
54.28+0.81 



47.62+2.36 
33.6113.29 
48.5912.49 



44.43+0.67 
46.77+0.78 
54.3110.81 



4.23 
7.23 
-1.91 



2.32 
4.27 
0.02 



7.39 



<0.01 
.08 <0.001 
5.3 n.s. 



12.96 
9.72 
9.47 



<0.05 
<0.01 
n.s. 



No 



* T-statistics has been estimated with a Welch's test using unequal variance. 

# Number of known SNPs in PCR product. 
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Figure 4 Verification of differential methylation of four arbitrarily chosen genes, a. Methylation differences between Red Junglefowl and 
White Leghorn in three different tissues, estimated by methylation sensitive high resolution melting (MS-HRM) analysis. The table shows the 
normalized fluorescent intensities together with the statistical analysis of 6-8 individual samples per population, at a temperature where positive 
and negative methylation controls showed highest intensity difference (T md ). The right column shows the number of SNPs present in the PCR- 
product amplified by the primers for the region analysed, b-d, The average methylation level of the promoter regions of three of the genes, 
FUCA1 (b), PCDHAC1 (c) and RUFY3 (d), in blood samples from F8-generation intercross birds between Red Junglefowl and White Leghorns. The 
birds differed as indicated in their genotype on an SNP-marker close to the differentially methylated locus. 



gene ontology (GO) analysis. We analysed the DM and 
the DE genes in each generation separately, and then 
selected those GO-terms and KEGG pathways (P < 0.1), 
which were significantly enriched in both generations 
(Additional file 6). 

A majority of the enriched GO terms are related to 
phosphorylation and kinase activity, important aspects 
of intercellular signalling. Looking specifically at the 
KEGG pathways enriched among DM and DE genes in 
offspring only (where the biological sample is consider- 
ably larger), the analysis shows that MAPK signalling 
pathway (which, for example, is associated with stress 
responses), long-term potentiation (affecting memory 



consolidation), neurotrophin signalling pathway 
(involved in neural differentiation) and GnRH signalling 
pathways (related to reproduction) are enriched. All 
these are potentially interesting from a domestication 
perspective, in that they may be related to well docu- 
mented differences between RJF and WL in stress toler- 
ance, behaviour and reproduction. 

Over-representation of epigenetically affected genes in 
selective sweep regions 

We considered that the epigenetic differences between 
the layer breed and their ancestor could reflect general 
effects of selection during domestication, as suggested 
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above, perhaps being related to differences in the 
domestication induced phenotypes, such as growth, 
feeding behaviour and social tolerance. If so, we would 
expect the epigenetic differences to be accumulated in 
genomic regions which have been under selection dur- 
ing domestication. Therefore, we compared our data to 
one of our earlier, and recently published, datasets on 
chickens [18]. This dataset consists of an extensive list 
of selective sweeps related to chicken domestication, 
based on resequencing of populations of RJF and a 
number of domesticated breeds. In total, 149 selective 
sweeps present in all domestic chickens, and 134 pre- 
sent in egg laying breeds only, were used. A sweep was 
defined as a 40 kb region where the heterozygosity 
Z-score was below -4. 

There were 216 DE genes (DE in both generations) 
with annotated loci within the 975 Mb of the genome 
covered by the sweep analysis. Five of them were situ- 
ated within 50 kb of selective sweeps present in all 
domestic chickens (non-significant association, based on 
a permutation test; P > 0.1), and nine in the laying 
breed sweeps (significantly more than expected by 
chance; P < 0.05). 

We performed the same analysis on 134 DM loci, and 
found that four were within 50 kb of sweeps in all 
domestic chickens (non-significant; P > 0.1), and six in 
laying breed sweeps (P < 0.05). The significant overlap- 
ping genes in laying breed sweeps are shown in Addi- 
tional file 7. 

It is interesting to note that ABHD7, which was the 
strongest DM and one of the strongest DE genes in our 
experiment, is positioned in a laying breed sweep. This 
gene is named EPHX4 in humans, and is related to 
detoxification of exogenous chemicals [19]. Based on its 
position in a selective sweep, and its differential methy- 
lation and expression, it would appear that the epige- 
netic variant of the gene (or the genotype affecting the 
epigenetic state of it) may have been selected during 
domestication. KSR1, an important gene in MAPK/Ras 
dependent signalling [20], as well as ADRA2C, an alpha 
adrenoreceptor that may be related to egg laying [21] 
and regulation of the sympathetic stress reaction [22], 
are also situated in laying breed sweeps. 

Although our data do not allow us to conclude on 
which genes and which sweeps that are associated with 
specific phenotypes, they suggest that selection of epige- 
netic variation may have been an important part of 
chicken domestication. 

General discussion 

Our findings show that differential methylation and gene 
expression in hypothalamus/thalamus are abundant in a 
comparison between domesticated White Leghorn layers 
and their wild ancestors, the Red Junglefowl. Many of 



these epigenetic differences are inherited, demonstrating 
transgenerational stability. It is possible that these differ- 
ences are a result of selection during domestication, tar- 
geting either sequence differences which affect epigenetic 
states of specific loci, or epigenetic states which are not 
related to sequence differences. 

The causal relationship between methylation and gene 
regulation is not clear, since differential methylation was 
associated with both up- and down- regulation of the gene 
expression, or did not affect it at all. Since similar dissocia- 
tion between methylation and gene expression has 
recently been found in the human genome as well [17,23], 
this indicates that epigenetic regulation is more complex 
than previously assumed. Whereas it is often believed that 
methylation of promoter regions is associated with down- 
regulation of gene expression, our results indicate that 
gene regulation is more complex than so. For example, 
chromatin structure may be more important than com- 
monly assumed. Furthermore, we found that CpG-islands 
are not always methylated, so there may also be evolution- 
ary contraints on methylation sites, hence affecting the 
rate with which epigenetic adaptations may occur in dif- 
ferent parts of the genome. Although speculative, these 
issues should be considered in future research. 

Some of the methylation differences observed appear to 
be tissue-specific, whereas others affect a wider range of 
cells. The mechanism whereby differential methylation at 
a particular locus only in, for example, the brain can be 
transferred from parents to offspring remains elusive. In 
Drosophila, similar observations have been made with 
respect to gene expression, where induced differences spe- 
cifically in the brain can be transmitted via sperm and 
cause tissue specific effects in the next generation [24]. 
Possibly, microRNA regulation may be involved [25], and 
both sperm and eggs may also transfer specific histone 
variants [26]. There is also a close link between genetics 
and epigenetics, in that the epigenetic state of a particular 
locus is determined by both genetic and epigenetic varia- 
tions at other loci [23,25]. 

Stable inheritance of epigenetic variants has been 
demonstrated in plants [7,13,27,28]. Also in mammals 
(mainly in rodents and humans), there is increasing evi- 
dence that this occurs widely [10,29]. Our results are the 
first to demonstrate the same in birds, and furthermore 
show a long-term stability over several generations of 
specific methylation states. 

Although we have only studied one population each of 
Red Junglefowl and domesticated chickens, the observa- 
tions in this experiment could indicate that selection of 
favourable epigenomes, or genotypes favourably affecting 
the epigenome, may have been an important aspect of 
chicken domestication. However, further studies are 
needed, where methylation of specific genes are analysed 
in a wide range of domesticated populations, analogous 
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to the recent study of sequence variation in the domesti- 
cated chicken genome [18]. 

Conclusions 

Both gene expression and promoter methylation profiles 
in hypothalamus differed between White Leghorn chick- 
ens and their wild ancestors, the Red Junglefowl Family 
differences within breed as well as breed differences 
were maintained in the offspring, showing a reliable 
inheritance of epigenetic states, which may have been an 
important mechanism involved in the rapid evolutionary 
changes of chickens during domestication. 

Methods 

Ethics statement 

The experiment and all its procedures were approved by 
the regional Ethical Committee. 

Animals and selection 

We used two breeds of chicken: a domesticated egg layer 
White Leghorn (WL) and a wild type Red junglefowl 
(RJF). The population backgrounds have been described 
elsewhere [30]. The WL was an outbred mixture of differ- 
ent leghorn breeds, established 1970, and since then kept 
in a closed population at the university. The RJF stem 
from an outbred zoo population, kept at the university for 
10 years (one generation per year) with a generation size 
of about 100 individuals. For more details about breeding 
and housing routines, see [6]. 

From each breed, two breeding pairs were selected based 
on divergent performance in four different fear tests con- 
ducted between 153-168 days of age. The selection criter- 
ion in each test was the frequency of standing/sitting alert, 
where a longer duration signifies a higher fear level [6]. 
The tests were: 1) Behaviour in an open field during 10 
min (this test was repeated at 286 days for check of consis- 
tency); 2) Behavioural response for 5 min after novel 
object introduction; 3) Behaviour 5 min before and 5 min 
after exposure to an aerial predator model; 4) Fearful 
behaviour toward a human. All tests were conducted in 
the same arena (0.5 x 1.5 m), and details of the test proce- 
dures have been described elsewhere [6]. All birds were 
weighed at 373 days of age. 

Eggs were collected for five weeks from each individual 
female in the breeding pairs and incubated in the same 
incubator in three consecutive batches. Numbers of eggs 
in the batches were balanced for family. Each batch con- 
sisted of 9-18 birds, and the batches were housed sepa- 
rately in mixed groups under the same conditions as the 
parents. At day 20 all offspring were tested in a similar 
open field arena as the parents. 

For HRM verification of the breed differences in methy- 
lation levels, we bred 8 offspring of each breed (four 
females and four males), using different animals (but same 



populations) as those above. These chicks were culled and 
sampled in the same way as described below. 

For studying the long-term stability of methylated loci, 
we used 18 birds from generation eight of an intercross 
between RJF and WL. The details of this intercross has 
been described elsewhere [3]. Briefly, one male RJF and 
three WL (same populations as those used in other parts 
of the present experiment) were used to breed 36 Fl, and 
these were intercrossed to produce about 1000 F2. From 
generation F3 onwards, about 100 birds per generation 
have been maintained by random mating and pedigree 
hatching up to generation F8. 

RNA and DNA isolation 

Parents were culled at day 373 after hatch and their off- 
spring at day 21. A sample of six male and six female off- 
spring (balanced between batches) were chosen from 
each family. A part of the brain enriched of thalamus/ 
hypothalamus was anatomically dissected and immedi- 
ately snap frozen in liquid nitrogen, and blood was 
collected immediately after culling. Samples were homo- 
genized in TRI-reagents (Ambion) using the FastPrep® 
-24 homogenization system with Lyzing matrix D tubes 
(MP Biomedicals). 

RNA was further extracted with the same method as 
has been used previously [8,9] following the protocol of 
the TRI-reagent manufacturer, except for a modification 
adding 0.25 ml isopropanol and 0.25 ml RNA-precipita- 
tion solution (1.2 M NaCl, 0.8 M disodium citrate). 
After the TRI-based RNA extraction, an RNeasy kit 
(Qiagen) was used to further purify the samples. 

DNA was extracted from the same TRI homogenate 
as the RNA. Precipitation was done on ice by adding 
150 \A of 100% ethanol to 300 \i\ TRI homogenate fol- 
lowed by gentle vortexing and 5 min 4°C incubation. 
After centrifugation in 12000x g for 10 min, the super- 
natant was discarded and the pellet resuspended in the 
ATL buffer of the DNeasy kit (Qiagen). 4 [i\ of 10 mg/ 
ml Rnase A was then added, incubated 2 min at room 
temperature, followed by addition of 20 |il of 10 mg/ml 
proteinase K (Qiagen) and 3 min incubation at 56°C. 
DNA was then extracted according to the DNeasy pro- 
tocol for animal tissue. Quality and quantity of both 
RNA and DNA was measured on a Bioanalyzer® instru- 
ment (Agilent Technologies) and NanoDrop® ND-1000 
spectrophotometer (Thermo Scientific). 

Both RNA and DNA extractions were treated indivi- 
dually for the parents, but as pools of six same-sex sam- 
ples in the offspring. 

Gene expression microarray 

This part of the experiment was performed at Uppsala 
Array Platform at Uppsala University, Sweden http:// 
www.medsci.uu.se/klinfarm/arrayplatform. A total of 16 
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GeneChip Chicken Genome Arrays (Affymetrix Inc.) 
were used to measure the expression of 33457 tran- 
scripts. Biotinylated fragmented RNA was prepared for 
each sample using standard procedures in GeneChip ® 
3' IVT Express Kit Users manual (Affymetrix Inc., Rev. 
1, 2008). This was followed by array hybridization for 
16 h in 45°C under constant rotation. Washing and 
staining was performed in a Fluidics Station 450 
and scanned using the GeneChip Scanner 3000 7 G 
(Affymetrix Inc.). 

Gene expression data analysis 

Analysis of the gene expression was performed using the 
statistical software R http://www.r-project.org with 
Bioconductor packages http://www.bioconductor.org. 
Normalization was done with the RMA method [31] 
and differentially expressed genes were evaluated using 
fold change in combination with a Bayes moderated 
t-test [32] adjusted for false discovery rate [33]. 

Correlation analysis, comparing the fold change of dif- 
ferential expression in parents and offspring, was done 
using Statistica v 9.1. Cluster analysis and heat maps 
were performed with the Genesis software v 1.7.5 [34]. 

DNA-methylation tiling array design 

For methylation analysis, we selected the genes which 
had the highest fold change of expression in both gen- 
erations of the breed comparison in the microarray data. 
From each gene, the promoter regions, defined as 7.25 
kb upstream and 3.25 kb downstream of the transcrip- 
tion start site (Ensemble genebuild WASHUC2), were 
used to design a custom 385 K DNA-methylation tiling 
array (Roche-NimbleGen). In total 3623 promoter 
regions were tiled to the array, with 50-75 mer probes 
and 100 bp median spacing, by the Madison design 
team at Roche-NimbleGen. 

Methylated DNA Immunoprecipitation (MeDIP) 

Protocols of MeDIP with buffer descriptions and general 
procedures have been published elsewhere [35,36]. Frag- 
mentation of 6 (ig thalamus/hypothalamus DNA was per- 
formed using a BRANSON sonifier 250 with a 13 mm 
disruptor horn (101-147-037) and a 3 mm tapered micro- 
tip (101-148-062). Samples were diluted with 450 ul 1 x 
TE in 1.5 ml tubes and sonicated at 10% amplitude by 
short 0.5 sec pulses (20 in total) with a rest between pulses 
of 0.5 sec. Fragment lengths of between 300-1000 bp were 
verified on a Bioanalyzer (Agilent technologies). After 
sample denaturation 10 min at 95°C, reference samples of 
10 [i\ was taken from each of the original samples and fro- 
zen. The remaining samples underwent methylated DNA 
immunoprecipitation by first diluting them in 1 x TE to 
450 ul, adding 51 ul of 10 x IP buffer and 10 ug of 5-meC 
antibody (Diagenode). Samples were then incubated in 



4°C for 2 h on a rotating platform, whereby 50 ul of clean 
Dynabeads Protein G (Invitrogen) in 1 x IP buffer was 
added and followed by an identical 2 h incubation. The 
beads-antibody-antigen complex was washed 3 times by 
placing the samples on a DynaMag-spin magnet, discard- 
ing the supernatant and adding 1 ml of 1 x IP. Complex 
digestion was done by adding 250 ul of Proteinase K diges- 
tion buffer and 5 ul Proteinase K (20 mg/ml), followed by 
rotation over night in 50°C. DNA was further purified by 
phenol/chloroform procedures with glycogen/ethanol 
-80°C precipitation. The pellet were washed in 100% etha- 
nol and resuspended in 60 |Lil 1 x TE. All samples, refer- 
ences as well as MeDIP's, were then whole genome 
amplified using the WGA2 kit (Sigma-aldrich) and puri- 
fied with QIAquick PCR purification kit (Qiagen). 

DNA-methylation tiling array labeling and hybridization 

Labeling and hybridization was performed at Roche-Nim- 
bleGen service lab at Iceland, Reykjavik, using standard 
protocols [37]. In short, MeDIP and reference samples 
were labeled with Cy5 and Cy3 respectively, using the 
NimbleGen Dual-Color DNA Labeling Kit. The MeDIP- 
Cy5 and reference-Cy3 samples from each tissue sample 
were then co-hybridised to the DNA-methylation tiling 
array using the NimbleGen Hybridization Kit and Hybridi- 
zation System. After washing with NimbleGen Wash Buf- 
fer Kit the slides were scanned by a NimbleGen MS 200 
Microarray Scanner. 

Tiling array data analysis 

Methylation data analysis was performed using Biocon- 
ductor in the open source R statistical software environ- 
ment [38]. To not loose genome wide methylation 
differences, the RINGO package [39] was used to prepro- 
cess the data within arrays by Tukeys biweight normaliza- 
tion and between arrays with A-quantile normalization. 
Significantly differentially methylated probes (FDR 
adjusted P-values) were extracted using the limma package 
[32]). Since promoters sometimes involved more than one 
significant probe, in all comparisons with the gene expres- 
sion data and correlations across generations, only the 
most significant probe of each promoter was considered. 
All significant probes that were stable across generations 
were checked for the occurrence of SNPs using a list of 
SNPs detected in a multibreed resequencing study recently 
published 18 

Methylation sensitive high resolution melting analysis 

Verification of differentially methylated genes, and analysis 
of differential methylation in alternative tissues and in 
F8-intercross birds, was done by methylation sensitivehigh 
resolutions melting (MS-HRM) analysis, principally 
as described by [40]. If not said otherwise, all procedures 
followed manufactures recommendations. DNA was 
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prepared from brain tissues as above and from blood using 
the DNeasy Blood and Tissue kit adjusted for nucleated 
blood (Qiagen). 

Positive control samples were synthesized by in vitro 
methylation, using a nuclease-free water diluted reaction 
mix of 16.5 (il, including an all bird pool of 1 (ig DNA, 
2 [i\ lOx NEBuffer2, 2 \i\ SAM (640 |iM), 1 [i\ SssI 
methylase (4 U/\il) (New England BioLabs Inc.). After 2 
h of incubation at 37°C, an additional 2.5 \i\ SAM was 
added to each sample, followed by another 2 h incuba- 
tion and then termination by heating at 65°C for 20 
min. 

Negative control samples were synthesized by whole 
genome amplification on the same all bird DNA pool 
(10-20 ng/(il) as for the positive control using the 
REPLI-g Mini Kit (Qiagen). The whole volume of ampli- 
fied negative controls were then mixed with 200 [il Buf- 
fer AL and 200 \i\ ethanol (99%) and purified with the 
DNeasy Blood and Tissue kit (Qiagen). 1 \ig DNA from 
both individual samples and controls were bisulfite trea- 
ted using the EpiTect Bisulfite Kit (Qiagen). 

PCR and High resolution melting was performed on a 
Rotor-Gene 6000 thermocycler (Corbett Research). 1 [i\ 
of the bisulfite treated samples/controls were prepared 
in a 10 [A PCR mix using EpiTect HRM PCR Kit (Qia- 
gen). A calibration series was also amplified using a 
mixture of positive and negative controls at 100%, 75%, 
50%, 25% and 0% of methylated DNA. PCR was per- 
formed in 45 cycles as follows: denaturation 10 s at 95° 
C, annealing at 30 s 54-55°C (primer dependent) and 
extension 10 s at 72°C. MS-HRM was run in the interval 
of 70°C to 90°C, with a 2 s 0.05°C steps, acquiring fluor- 
escence data at the Rotor-Gene HRM channel. Primer 
sequences and annealing temperatures can be seen in 
Additional file 8. All MS-HRM reactions were run in 
triplicates. 

Annotation and GO analysis 

Affy Chicken ID, EntrezGene ID, EnsembleGene ID, 
WikiGene ID and chromosomal regions were extracted 
from the Affymetrix annotation file (release 29), and 
further annotated with Ensemble's BioMart tool [41]. 
We used DAVID 6.7 http://david.abcc.ncifcrf.gov to 
extract significantly enriched gene ontology terms and 
KEGG pathways within our datasets [42,43]. To increase 
the possible DAVID hits we first extracted the homolo- 
gous human Ensembl ID of our chicken genes in Bio- 
Mart [41]. CpG island prediction was performed with 
EMBOSS CpGPlot [44]. 

Analysis of sweep overlaps 

219 DE genes and 134 DM promoters (significant in 
both parents and offspring) fell within the 975 Mb that 
previously has been searched for selective sweeps 18 . To 



investigate whether genes or promoters were signifi- 
cantly enriched in sweep regions, 1000 sets of random 
intervals were generated over the 975 Mb for each ana- 
lysis, each interval in each set chosen to represent one 
DE or one DM gene. The number of overlaps between 
the randomly generated interval and a sweep (within 50 
kb of the sweep) was compared to the actual number of 
real overlaps. A probability of the observed coincidence 
of less than 5% was taken as a significant association 
between DM/DE genes and sweeps. 

Additional material 



Additional file 1: Phenotypes of parents and offspring. The 

behavioural scores and weight data for the animals in the experiment. 

Additional file 2: Nrs of differentially expressed and methylated 

genes. Total numbers of significantly differentially expressed genes or 
methylated promoters (P < 0.05, FDR-corrected) in the different 
comparisons. 

Additional file 3: Differentially expressed genes. A full list of all genes 
found to be differentially expressed, comparing breeds, in both 
generations, including their chromosomal alignment and accession 
numbers. 

Additional file 4: Expression and methylation. Gene expression 
differences plotted against promoter methylation differences between 
WL and RJF offspring. 

Additional file 5: Differentially methylated genes. A full list of all 
genes, where the promotors were found to be differentially methylated, 
in both generations. 

Additional file 6: Gene function. Gene ontology and KEGG pathway 
annotation for the genes which were either differentially expressed or 
differentially methylated in both generations, comparing between 
breeds. 

Additional file 7: Selective sweep representation. Genes which were 
differentially expressed or methylated between breeds in both 
generations, and significantly overrepresented in selective sweeps 
associated with domestication. 

Additional file 8: Primer structures. The bisulfate converted primers 
used in HRM. 
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